Influence of Material Parameter Variability on the Predicted Coronary Artery Biomechanical Environment via Uncertainty Quantification

Central to the clinical adoption of patient-specific modeling strategies is demonstrating that simulation results are reliable and safe. Indeed, simulation frameworks must be robust to uncertainty in model input(s), and levels of confidence should accompany results. In this study, we applied a coupled uncertainty quantification-finite element (FE) framework to understand the impact of uncertainty in vascular material properties on variability in predicted stresses. Univariate probability distributions were fit to material parameters derived from layer-specific mechanical behavior testing of human coronary tissue. Parameters were assumed to be probabilistically independent, allowing for efficient parameter ensemble sampling. In an idealized coronary artery geometry, a forward FE model for each parameter ensemble was created to predict tissue stresses under physiologic loading. An emulator was constructed within the UncertainSCI software using polynomial chaos techniques, and statistics and sensitivities were directly computed. Results demonstrated that material parameter uncertainty propagates to variability in predicted stresses across the vessel wall, with the largest dispersions in stress within the adventitial layer. Variability in stress was most sensitive to uncertainties in the anisotropic component of the strain energy function. Moreover, unary and binary interactions within the adventitial layer were the main contributors to stress variance, and the leading factor in stress variability was uncertainty in the stress-like material parameter that describes the contribution of the embedded fibers to the overall artery stiffness. Results from a patient-specific coronary model confirmed many of these findings. Collectively, these data highlight the impact of material property variation on uncertainty in predicted artery stresses and present a pipeline to explore and characterize forward model uncertainty in computational biomechanics.


Introduction
Physics-based simulations of the cardiovascular system are increasingly being integrated into clinical decision-making (Douglas et al. 2015;Driessen et al. 2019), surgical planning (Trusty et al. 2019), and medical device design (Timmins et al. 2011).
Moreover, the U.S. Food and Drug Administration (FDA) has published widely on using simulations to promote the safety, effectiveness, and security of FDA-regulated products (Morrison et al. 2017(Morrison et al. , 2018;;Pathmanathan et al. 2017;Ahmed et al. 2023;Food and Drug Administration: Center for Devices and Radiological Health 2023).As simulations contribute to clinical workflow and regulatory approval and may affect downstream outcomes (e.g., major adverse events, patient death), there is a pressing need to provide confidence in simulation predictions and demonstrate that results are reliable and safe before clinical adoption.Such confidence in simulation pipelines is available in idealized scenarios but is marred by uncertainties, which manifest through variability in the subject and clinical variabilities (i.e., model input parameters) that cloud the predictive and prognostic lenses of computer-based modeling.Central to the clinical adoption of patient-specific modeling strategies, therefore, is clearly demonstrating that simulation results are reliable and safe.That is to say, it is essential that simulation results be accompanied by levels of confidence when they potentially impact life-altering decisions.
Advances in medical imaging, computational mechanics, biomechanics, and computing power now enable simulations that predict arterial tissue deformations at the patient-specific level (Taylor and Figueroa 2009;Taylor and Humphrey 2009).In addition to the model geometry, boundary conditions, and numerical approaches to solve the governing equations, the constitutive relation(s) describing the behavior of the material(s) under conditions of interest are required to compute the transmural wall stresses that influence the homeostatic and maladaptive mechanobiological processes.
While experimental approaches have been developed to characterize the non-linear, pseudoelastic, and anisotropic material response of vascular tissue under loading, there is much variability in the employed techniques.For example, such methods as ring tests, in-plane biaxial tests, pressure-diameter tests at the in vivo length, and biaxial tests consisting of cyclic pressure-diameter and axial force-length protocols have been utilized to characterize the mechanical properties of vascular tissue (see the comprehensive review by Feruzzi et al. and references within (Ferruzzi et al. 2013)).
Regardless of the specific form of the strain energy function (SEF, !), regression analysis can identify the best-fit values of the material parameters within.In addition, there is variability in the applied regression method (e.g., Marquardt-Levenberg) and candidate objective function that is minimized (Humphrey 2002;Ferruzzi et al. 2013).
Due to the variability within experimental testing protocols and fitting approaches, as well as variability within and across tissue samples, there exist inherent uncertainties in material parameter(s) describing the soft biologic tissue that propagate to the simulation-predicted mechanical environment.
In the present study, we incorporated advancements in the field of UQ to evaluate the variability in the output of computational simulations of the arterial mechanical environment due to intrinsic uncertainty in material parameter estimation.In contrast to a traditional deterministic simulation where input parameters have a fixed value that results in a single model output, UQ provides a statistically rigorous approach to determine the influence of input parameter uncertainty by examining a distribution of model outputs (Najm 2009).We applied a novel open-source UQ software tool, UncertainSCI, which employs polynomial chaos expansion (PCE) to assess sensitivity, to a forward-modeling framework (Narayan et al. 2022).Therefore, the goal of this study was to leverage PCE UQ to examine the impact of uncertainty in tissue material properties on the variability in model outputs, namely the predicted stress under physiology loading.Given the clinical significance of coronary artery disease and the role of mechanics in the development and progression of the disease (Brown et al. 2016;Tsao et al. 2022), we focused on uncertainty in material characterization and computational models of this vascular territory.To demonstrate the approach, we evaluated models of a generalized multi-layered, thick-walled vessel representative of a coronary artery and a patient-specific model of an epicardial coronary artery.

Methods
An overview of the integrated UQ-finite element (FE) modeling framework is presented in Fig. 1.Briefly, probability distributions were fit to "-material parameters, which were derived from material testing of human coronary tissue, in a structurallymotivated SEF.The parameter space was sampled to generate #-parameter ensembles.Utilizing a batch-processing framework, forward FE models for each parameter ensemble were created, and FE analysis was carried out to predict tissue deformation, strains, and stresses.Finally, statistics of the model outputs were computed, and uncertainty was quantified due to material parameter variability.

Material Parameter Probability Distribution Sampling
Material parameters derived from previous layer-specific mechanical behavior testing of 13 human nonatherosclerotic left anterior descending coronary artery tissues were employed herein (Holzapfel 05).Briefly, uniaxial extension testing was performed in the circumferential and longitudinal (axial) directions for the media and adventitia layers.Best-fit material parameters were determined from the mechanical behavior curves using the structurally-motivated SEF, where 4 represents the ground matrix stiffness, 0 1 is a fiber stress-like parameter, 0 2 is a dimensionless parameter, 5 is a measure of fiber dispersion within the bounds of [0,1] (0 = no fiber alignment, 1 = perfect fiber alignment along the prescribed vector defined by angle 6), and + 1 and + 2 !are the first and fourth invariant, respectively, of the right Cauchy-Green tensor (%), defined as, for a deformation gradient that takes the form, < = =;>?(7 3 , 7 4 , 7 5 ).' ! is a unit vector, $0, 89:6 !, :;"6 !(, indicating the orientation of a fiber family, where 6 !defines the angle between the embedded fiber family and the circumferential axis in the circumferentialaxial plane.Two fiber families were considered (6 " = −6 # ), together describing symmetric fiber families with the same material properties around the long axis of the vessel.
Probability distributions were created for each constitutive parameter in the medial and adventitial layers (10 in total, Fig. 2).Note that parameter distributions were assumed independent, which was guided by the original description of the SEF (Eq. 1) and a lack of data demonstrating any physical relationship amongst them (Holzapfel et al. 2000).Material parameters 4, 0 1 , 0 2 , and 6 employed gamma distributions to avoid non-physiologic parameters (i.e., values had to be >0), whereas 5 used a beta distribution to take advantage of its inherent bounds [0,1].Probability density functions (PDFs) were fit to the experimental data of each individual parameter using maximum likelihood estimation via the mle function in Matlab.The material parameter PDFs were defined as distribution objects within UncertainSCI, and a parameter ensemble was created from the 10 PDFs that sampled the entire parameter space and accounted for potential output dependence on parameter ensembles (Fig. 1a).
Figure 2 here.Fig. 2 Layer-specific material parameter distributions and probability density functions (PDFs).Material parameters from all samples for the media (top row) and adventitia (bottom row) layers from prior uniaxial testing of coronary arteries (Holzapfel et al. 2005).PDFs for each parameter were fit to the observed data.Gamma distributions were fit to data for parameters 4, 0 " , 0 # , and 6, and a beta distribution was fit to data for 5.

Idealized Artery Computational Model
A generalized computational model of a human left anterior descending coronary artery was constructed.The artery was modeled as a multi-layered, axisymmetric quarter-cylinder (A = 1 mm, B ! = 1.59 mm, B 6 = 2.25 mm), with a medial and adventitial layer thickness of 0.32 and 0.34 mm, respectively (Holzapfel et al. 2005).The intima layer was neglected in this model, as it provides negligible structural support (Burton 1954).The structurally-motivated SEF given by Equation 1, an available material model in FEBio (HGO-coronary), described both the medial and adventitial layers.The artery was discretized with 8-node hexahedral "brick" elements with 6 elements in the radial direction for each layer.A mesh convergence study demonstrated that 12 elements in the radial direction (6 elements in each layer) were required to achieve convergence, as a higher mesh density (24 radial elements) led to a <2.5% change in the 2-norm criteria for the 1 st principal stress (Supp.Fig. 1).The media and adventitia arterial layers were "welded" with shared nodes at the interface.Applied loads included lumen pressures of 80 (diastolic) and 120 (systolic) mmHg.The boundary conditions comprised of fixing the vessel ends in the axial direction and symmetry in the C-planes to restrict rigid body motions.Quasi-static finite element analysis was performed using the open-source, nonlinear finite element software suite FEBio (Maas et al. 2012(Maas et al. , 2017)).Solver details include using the implicit solver, auto-time stepper (initial and maximum time-step size = 0.1), non-symmetric form of the stiffness matrix, quasi-Newton method (Broyden-Fletcher-G-S; global stiffness matrix reformed each time step), and PARDISO linear solver.Solver settings ensured numerical robustness and the ability to support parallel execution.Simulation results were post-processed to evaluate the deformed inner (D ! ) and outer (D 6 ) radii, transmural distributions of 1 st and 3 rd principal stresses (E 1 , E 3 ), and distensibility (F, Eq. 2) (Ferruzzi et al. 2013), which was defined as, , (Eq. 2) where = !,9:9and = !,7!;9 are the deformed inner diameters at systole and diastole, respectively, and G 9:9 and G 7!;9 are systolic and diastolic pressure, respectively.
A batch-processing scheme was developed to iterate through material parameter ensembles while using the same FEBio input file, which contained information on model mesh and connectivity, material SEFs and parameters, and boundary conditions.
A Matlab sub-routine was written that iterated through the size-H parameter ensemble, writing H unique input files (Fig. 1b).FEBio was called and executed using the GIBBON toolbox (M Moerman 2018), and simulation results were stored for UQ analysis.For the idealized artery models, the batch-processing scheme and FE models were run on a Windows 10 server machine with an Intel® Xeon® Silver 4110 CPU (8 cores at 2.10 -3.00 GHz).

Patient-specific Coronary Artery Computational Model
A three-dimensional representative patient-specific model of the left main and left anterior descending coronary arteries was constructed by expanding established techniques (Samady et al. 2011;Timmins et al. 2015).The end-diastolic geometry was created by fusing bi-plane angiographic image data and virtual-histology intravascular ultrasound (VH-IVUS) images (Fig. 3).Lumen and media-adventitial boundary contours were stacked perpendicular to the IVUS catheter centerline, and catheter torsion was accounted for via the sequential triangular algorithm (Wahle et al. 1999).A medial layer was constructed from smoothed IVUS boundary contours, and an adventitial layer was added with a constant thickness of 400 μm (Waller et al. 1992).Branches were added from IVUS and angiographic-defined locations with branch layer thicknesses derived from post-mortem coronary mean lumen diameter and layer thickness values (Waller et al. 1992).The geometry was meshed with nonlinear tetrahedral elements via tetGen and Gibbon (Si 2015;Maas et al. 2016;M Moerman 2018), with unique material properties for each layer prescribed using structurally-motivated SEF (Eq. 1) (Holzapfel et al. 2005).To aid the application of boundary conditions, a rectangular box of perivascular (PV) tissue with compressible, neo-Hookean properties (I = 1 kPa, J = 0.3) was added around the coronary geometry and shared identical nodes with the outer vessel surface.The PV outer boundary surfaces were at least 10 mm away from all nodes in the artery (Fig. 3B) and were fixed in all global directions.Preliminary studies on idealized and patient-specific coronary geometries demonstrated that a PV support with those material properties (I = 1 kPa, J = 0.3) and a thickness of 10 mm had a negligible effect on the deformation of the arterial tissue under an applied lumen pressure.Axial motion was prohibited at the vessel and branch end surfaces.The lumen was pressurized to 40 mmHg (note: the reference geometry represented the vessel at end-diastole; ~80 mmHg).A set of patient-specific models was created using an identical set of material parameter ensembles from the ideal quarter cylinder model at PCE order 3 (n=592 parameter samples).Given the increased complexity of the patient-specific models compared to idealized models, solver and solution parameters were modified.Broyden's quasi-Newton method was employed as the solver.The autotime stepper parameters were reduced (initial time-step size = 0.01, max time-step size = 0.05), and the aggressiveness parameter was turned on to aid the identification of time-step size after a failed step.To aid computational efficiency and solution convergence, the discretization was refined until mesh quality, as determined by the radius-edge ratio in the Mesh Inspector feature in FEBio, was deemed suitable.Across the >169k quadratic tetrahedral elements in the patient-specific artery model, excluding perivascular support, the average radius-edge ratio was 0.934, and <0.003% of the elements had a value >2.5.The batch scheme was executed on two high-performance compute servers (192 Intel Xeon Platinum 8360H CPU @ 3.00GHz cores (HT) per machine) with 14 concurrent jobs running on each machine using 8 cores per simulation (Scientific Computing and Imaging Institute, University of Utah).

Uncertainty Quantification and Sensitivity Analysis
The open-source, Python-based software suite, UncertainSCI, was employed to perform forward model UQ analysis (Narayan et al. 2022).UncertainSCI utilizes non-intrusive PCE techniques to query forward model data output over a parameter ensemble (i.e., generated parameter samples) to construct a parameter-to-model-output emulator.With a tunable parameter K, the PCE order, the emulator is comprised of a sum of polynomial functions of degree at most K, and serves as a surrogate model that approximates the mapping between the input parameter(s) and model output without the need to solve for a computationally expensive forward model.In this study, for example, the constructed emulator provides a relationship between SEF parameters (inputs) and FE-predicted principal strains and stresses (outputs) across the queried range of parameter distributions.UncertainSCI constructs parameter samples via the weighted approximate Fekete points method (Guo et al. 2018).The construction of the emulator allows direct extraction of statistics, uncertainty characteristics, and model sensitivity.For a comprehensive description of polynomial chaos techniques and their application, the interested reader is directed to the work by Najm (Najm 2009).
The quality of the PCE was assessed in the idealized computational models by quantifying the relative error (L < ) between the PCE approximations and FE-predicted model outputs (e.g., E 1 ), whereby, where NO P approximates the solution to the model output (Q R⃑ ), T is the number of elements through the vessel thickness (i.e., transmurally), and U|⋅|U 2 indicates the 2norm of the vector.For orders K = {1,2, . . .,5}, the relative error was calculated for 5 independent PCE runs.In addition, parameter ensembles were oversampled to evaluate error stability across sampling rates and ensure the aliasing error is first-order Sobol indices) and parameter ensembles (binary, tertiary, etc. interactions) have on the variance in the model output for 1 st principal stress (E1).

PCE Construction, Quality, and Order Convergence
The number of parameter ensembles generated (i.e., H, Fig. 1b), time required to generate these ensembles within UncertainSCI, and run time for the FE-batch processing within FEBio for the idealized geometry across PCE orders are presented in Table 1.The number of parameter ensembles ranged from ~40 (order 1) to >6,000 (order 5) and required between several minutes to many hours to generate the ensembles and run the FE simulations.Evaluation of the relative error for the 1 st and 3 rd principal stresses across orders demonstrated reduced error and error variability across multiple PCE runs as order number increased (Fig. 4a,b).Data indicated that order 3 captured the dominant uncertainty modes, as relative error values were <0.4% with a standard deviation of <0.1 across 5 runs.Furthermore, Sobol indices maintained stability at order 3 and higher orders.First-order Sobol indices in the adventitia changed by <0.003 and the relative positions remained unchanged across orders 3 to 4 (Fig. 4c).
Smaller changes were observed at higher orders (Supp.Fig. 2).A similar trend in the stabilization of the first-order Sobol indices was seen in the medial layer (Supp.Fig. 2).
Moreover, second-order Sobol indices and their relative positions were preserved across orders 3 to 5 (Supp.Fig. 3).Examining the FE results from the order 3 simulations demonstrated that a range of deformed geometries and distensibilities were captured (Fig. 5).Across the ~600 parameter ensembles at order 3, deformed inner diameter and thickness values ranged from 3.76-4.78mm and 0.49-0.58mm, respectively, and distensibility values ranged from 4.67 to 18.72 MPa -1 .

UQ and Sensitivity Analysis in Idealized Coronary Model
The propagation of material parameter uncertainties to the transmural distribution of 1 st principal stress at 120 mmHg yielded large deviations from the median within the medial and adventitial layers (Fig. 6a).While median stress values and stress variance decreased radially through each layer, there was an abrupt increase in the variance at the innermost region of the adventitia.Also, variances were higher overall in the adventitia.As a result, coefficient of variation values in the adventitia were >1.5× the values in the media, indicating adventitial stress values had greater dispersion around the mean (Fig. 6b).Sensitivity analysis highlighted that the material parameters in the anisotropic component of the adventitia dominated the variance in the predicted 1 st principal stress (Fig. 7).For example, adventitial material parameters 0 1 , 0 2 , 6, and 5 accounted for nearly 70% of the variance in predicted stress values due to a single parameter (unary interactions), with 0 1 alone accounting for ~25% (i.e., variance in FEpredicted stresses are largely explained by the uncertainty in the stress-like parameter describing the contribution of the adventitial fibers to the artery stiffness, 0 1 ).Notably, the uncertainty in the stiffness of the Neo-Hookean ground matrix, controlled by parameter 4, had a negligible effect on stress variance in the media, and only a marginal effect in the adventitia.Sensitivity analysis further revealed unique interactions between two parameters (i.e., binary interactions) that contributed to the variance in predicted stress.While unary interactions were dominant, binary interactions still accounted for 12.2% of the variance in predicted 1 st principal stress (Fig. 8a).Examining pairwise interactions within each arterial layer highlighted that such interactions in the adventitia accounted for far greater stress variance than those in the media (Fig. 8b).Of the 12.2% of the variance in stress due to binary interactions, the interaction involving the adventitia alone accounted for 46.8% of the variance, and the media-only interaction accounted for 11.6%.The binary interaction of parameters across layers (i.e., inter-layer) accounted for 41.6% of that variance.That is to say, 5.7% of the (total) variance in predicted stress was due to binary interactions between adventitial parameters alone (46.8% of 12.2%), compared to 1.4% for medial parameters and 5.1% for inter-layer parameter combinations.Across the 45 possible pairwise combinations, interactions between the adventitia 6-5 and 0 # -6 dominated, with normalized second-order Sobol indices of 0.21 and 0.10, respectively (Fig. 8c).Moreover, 5 of the top 10 binary interactions were between adventitial SEF parameters within the anisotropic component.Modest interactions between parameters within the media and adventitia were observed (e.g., media 0 " , adventitia 0 # : 0.07), and only 1 binary interaction between medial parameters was in the top 10 (media 0 " , 0 # : 0.07).Lastly, tertiary interactions (i.e., interactions between 3 material parameters) accounted for <3.1% of the variance in predicted 1 st principal stress (Fig. 8a).

Application of UQ to Patient-specific Model of Coronary Artery
The batch-processing scheme (592 patient-specific models, PCE order 3) completed in ~49.5 hours on the multicore compute servers.Examining the uncertainty in 1 st principal stress revealed spatial heterogeneity in statistical and sensitivity measures across the physical domain (Fig. 9a).At a cross section distal from the left circumflex (slice 1; Fig. 9b), for example, mean stresses across the simulations ranged from 15.9-44.6 and 11.8-27.0kPa in the medial and adventitial layers, respectively.
Similar trends of higher values in the media were observed at other spatial locations (Fig. 9b) and when comparing standard deviation and coefficient of variation values across the models.Like the idealized model results, first-order Sobol indices in the patient-specific model associated with adventitial material parameters dominated variances in predicted stress (Fig. 9c).The isotropic and anisotropic material parameters in the adventitia accounted for 23.3% and 42.5% of the variance in stress, respectively.Notably, first-order Sobol indices have marked spatial variation, with increased dispersion in the adventitial layer (Fig. 9c).Given that vascular tissue has an anisotropic structural organization and a nonlinear stress-strain response, it was not surprising to see that uncertainty in the anisotropic SEF parameters accounted for the greatest variance in stress (Figs. 7,8).
Moreover, sensitivity analysis demonstrated that the unary interaction of the initial stiffness of the fibers (0 1 ) and the binary interaction between fiber angle (6) and fiber dispersion (5) in the adventitia layer were prominent parameters that influence stress variabilities.Experimental data highlight that the adventitia is stiffer than the media layer (Holzapfel 2005), resulting from the dense network of type I collagen within the ground matrix.Thus, the presented results provide further evidence of the influence of the arterial structure on the mechanical behavior of this soft tissue.Also, it is important to recognize the greater dispersion of the anisotropic SEF parameters in the adventitia versus media (Fig. 2), contributing to the larger coefficient of variation and Sobol indices in the adventitia.In addition to increased sample testing, advances in experimental approaches to better characterize and quantify these SEF parameters will promote reduced uncertainty in calculated stress values.
Efforts have been focused on integrating computational models of the cardiovascular system into the clinical setting; however, multiple sources of uncertainty must be accounted for to provide confidence in model predictions.In patient-specific models, for example, uncertainty arises in domain (geometry) construction, boundary conditions, numerical schemes, and, as investigated herein, material properties.In the context of material properties, sources of uncertainty are present when utilizing either population-based or patient-specific data.A population-based approach was used in this study, whereby variability across patient samples contributed to the uncertainty in arterial stiffness (Holzapfel et al. 2005).Importantly, this approach requires quantifying arterial mechanical properties on large data sets utilizing standardized protocols, including biaxial material characterization, to minimize experimental variance (Walsh et al. 2014).While methods to non-invasively determine patient-specific material properties are in development (e.g., DENSE MRI (Bracamonte et al. 2020)), there remain limitations with such strategies that must be resolved before adoption.Moreover, uncertainty exists in patient-specific strategies due to material spatial heterogeneity and image noise that must be accounted for within the modeling framework.Despite longitudinal studies that have utilized deterministic modeling approaches to demonstrate the utility of arterial wall stress as a prognostic marker of, for example, coronary plaque or abdominal aneurysm rupture (Teng et al. 2014;Polzer et al. 2020), effective modeling strategies and decision guidelines must be robust to uncertainty and provide levels of reliability and safety before clinical adoption.
Capturing the randomness of the input parameters (i.e., defining accurate PDFs) is central to parametric UQ, whereby the randomness of input parameters propagates forward to the model outputs.The accuracy of constructed PDFs to capture SEF material parameter distributions depends on whether the evaluated samples (observations) sufficiently explore the parameter space.While mechanical testing data on human tissue samples across the vascular tree are reported (Vande Geest et al. 2004;Holzapfel et al. 2005;Teng et al. 2009), these studies are often limited to a few samples, which may not provide sufficient numbers to define representative PDFs.This study fit PDFs to SEF parameters derived from the experimental testing of 13 human coronary tissue samples, which is the largest reported data set on layer-specific mechanical testing of human coronaries (Holzapfel et al. 2005).The fit PDFs capture the distribution of observations (Fig. 2); however, it is unclear if these fits represent the population distribution and whether the selected distributions are the best descriptors.
Indeed, a study reported biaxial testing data from 125 human femoropopliteal arteries and provided a comprehensive analysis that examined differences in material properties across age and disease severity (Jadidi et al. 2021).Yet, given the difficulty in procuring human tissues, particularly healthy samples, such studies are rare, especially in the case of coronary arteries.Thus, the standardization of testing protocols is warranted to allow consolidation of data sets towards improved characterization of population distributions and quantification of parametric UQ input parameter randomness.
An advantage of UncertainSCI is that it utilizes non-intrusive UQ techniques, which do not require changes to existing simulation frameworks or numerical schemes, to calculate accurate statistical measures of the forward propagation of uncertainty.
Moreover, UncertainSCI uses PCE techniques, which are more efficient and offer better convergence than Monte-Carlo (MC) and quasi-MC-based approaches (Xiu and Hesthaven 2005) and offer advantages in biomedical simulations, where the dependence on parameters is often smooth.Although PCE and MC approaches sample the multivariate parameter distributions, MC approaches require more samples to obtain reliable sensitivity measures and are thus more computationally demanding (Eck et al. 2016).Regardless of the sampling approach, and particularly relevant to soft biological tissue, it must be ensured that the sampled parameter ensembles produce physical stress-stretch responses (Robertson and Cook 2014).While MC methods have been successfully applied to cardiovascular simulations (Sankaran and Marsden 2011;Tran et al. 2019), these approaches can be problematic for complex patient-specific models (e.g., Figs. 3 and 9).Alternatively, PCE approaches compute equivalent statistical metrics, with orders of magnitude fewer evaluation samples compared to MC approaches (Eck et al. 2016;Burk et al. 2020).Importantly, however, PCE approaches are only recommended when the number of uncertain parameters is limited, typically less than 20, after which PCE strategies are no longer more efficient than MC methods (Xiu and Hesthaven 2005;Crestaux et al. 2009;Eck et al. 2016).While the presented study evaluated 10 parameters describing the material properties for the UQ analysis (Eq.1, Fig. 2), additional areas of uncertainty are present in cardiovascular simulations, as discussed above, that would increase UQ complexity and computational demand.
Methods can be employed to reduce the number of uncertain inputs.For example, if uncertain inputs have minimal effects on model output variance (i.e., small first-order Sobol indices, Fig. 7), those inputs can be fixed within their uncertain domain.Recently, a novel UQ framework that utilized a multilevel multi-fidelity MC estimator, which incorporates results from zero and one-dimensional models across mesh (spatial) resolutions to efficiently construct estimators, was shown to greatly reduce computational costs (10 to 100× reduction) for UQ in hemodynamic simulations (Fleeter et al. 2020).Continued advancements in data-efficient UQ methods to promote clinical translation and adoption are warranted.
The lack patient-specific material properties brings into question the reliability of existing coronary artery model results.Without directly comparing results derived from patient-specific or population-based material properties, a few remarks can be made regarding reliability.First, advances in medical imaging, hemodynamic assessment, and constitutive modeling have promoted patient-specific modeling capabilities beyond 2D, linear elastic computational models derived from histology data (Cheng et al. 1993).Thus, even with the lack of material property data, patient-specific models better represent the in vivo anatomy and physiology, and model results are more reliable.Second, even with material property limitations, current modeling efforts enable hypothesis generation and testing that have revealed new insights into cardiovascular biology and medicine.As an example, early observations and hypotheses on the role of plaque stress in coronary plaque rupture motivated studies that have realized these observations in patients with acute coronary syndrome (Richardson et al. 1989;Teng et al. 2014).Third, the presented UQ analysis indeed provides an assessment of the reliability of modeling results given material property uncertainty.Thus, studies can (and should) address the interpretation of their results within the context of our presented findings.What remains unknown is how the uncertainty in these data translates to correlations with clinical observations or outcomes, promoting the predictive power of biomechanical indices.
The presented results have implications for interpreting correlations between biological processes, such as tissue homeostasis, growth and remodeling, and disease progression, and the mechanical environment, such as stresses and strains (i.e., vascular mechanobiology) (Humphrey and Schwartz 2021).Whether utilizing analytical solutions or FE approaches to determine stresses within thick-walled vascular tissue, variability in calculated/predicted stresses due to material parameter uncertainty could impact the drawn correlations.Moreover, stress mediated growth laws (i.e., constitutive relations that describe cellular and extracellular matrix produce and removal rates as a function of stress) are central to constrained mixture models of soft tissue growth and remodeling (Humphrey 2021).Recognizing and quantifying the impact of uncertainty on these relationships is critical to advance the understanding of the evolution of soft tissue geometry, composition, and material behavior under complex loading.
There are limitations in this study that should be acknowledged; however, these limitations do not detract from the significance of the results.First, not all possible PDFs were explored to describe the distribution of the reported SEF material parameters.
While PDFs that yielded non-physical material parameters were excluded, continued investigation of PDFs that best describe the experimental data is warranted.Second, residual strain was not included in the computational models.While residual strain homogenizes the stress field in the artery wall under physiologic loading (Chuong and Fung 1986), these do not mitigate the variance in wall stresses due to uncertainty in material parameters.Importantly, the opening angle, the geometric quantity that aids quantification of the displacement field to determine residual strain, is yet another variable with uncertainty to evaluate and determine its impact on arterial stress field variance.Third, material testing data on the intimal layer reported in the Holzapfel et al.
study (Holzapfel et al. 2005) was not incorporated into the idealized model.These data were derived from the mechanical testing of specimens with non-atherosclerotic intimal thickening and diffuse intimal hyperplasia, which are detectable with VH-IVUS imaging (García-García et al. 2009).Future investigations utilizing the presented FE-UQ framework will seek to incorporate additional tissue components and plaque phenotypes.Lastly, only one patient-specific geometry was explored.Although the coronary anatomy and spatial variation in sensitivity measures impart complexities that make it difficult to draw immediate conclusions, the demonstrated application of the UQ-FE framework provides a novel approach for future investigations.

Conclusion
In summary, we present a computational framework to explore, characterize, and quantify forward model uncertainty in FE simulations of the arterial wall biomechanical environment.We report that uncertainties in SEF parameters describing the material response of a multi-layered, thick-walled artery under physiologic loading are pushed forward, leading to considerable variances in transmural stress fields.These data highlight that there remains a pressing need to promote experimental data collection towards better characterizing SEF material parameter distributions, and further understanding the propagation of such uncertainty to the predicted kinematics and stresses.Moreover, our efforts demonstrate the demand for continued rigor in computational biomechanics by providing confidence in calculated stress metrics to address complex biological and clinical problems.Fig. 2 Layer-specific material parameter distributions and probability density functions (PDFs).Material parameters from all samples for the media (top row) and adventitia (bottom row) layers from prior uniaxial testing of coronary arteries (Holzapfel et al. 2005).PDFs for each parameter were fit to the observed data.Gamma distributions were fit to data for parameters 4, 0 1 , 0 2 , and 6, and a beta distribution was fit to data for 5.       as mean ± standard deviation across all elements within each material layer.

Figure
Figure 3 here.

Fig. 3
Fig. 3 Patient-specific coronary model construction.a Angiographic and VH-IVUS data

minimized 1 .
Statistical measures (e.g., mean, standard deviation, coefficient of variation) were calculated directly from the PCE model output.Sensitivity indices, which measure the relative contribution of individual parameters and parameter ensembles to the overall variability of the emulator (i.e., variability in model output), were calculated across the SEF parameters.More specifically, Sobol indices (Sobol′ 2001) were determined to measure the direct effect of an individual parameter (unary interaction;

Fig. 4
Fig. 4 PCE and sensitivity indices convergence.Relative error (L < ) in a 1 st principal stress and b 3 rd principle stress for PCE analysis across orders and runs.c First-order

Figure
Figure 8 here.

Fig. 8
Fig. 8 Second-order sensitivity analysis from order 3 PCE analysis.a Percent output

Figure Captions Fig. 1
Figure Captions

Fig. 3
Fig. 3 Patient-specific coronary model construction.a Angiographic and VH-IVUS data

Fig. 4
Fig. 4 PCE and sensitivity indices convergence.Relative error (L < ) in a 1 st principal stress and b 3 rd principle stress for PCE analysis across orders and runs.c First-order

Fig. 5
Fig. 5 Distribution of deformed geometries and structural stiffness from FE models for

Fig. 8
Fig.8Second-order sensitivity analysis from order 3 PCE analysis.a Percent output

Fig. 9
Fig. 9 Uncertainty quantification and sensitivity analysis in a patient-specific coronary

Table 1
Computational times across PCE orders for the idealized artery model.Five PCE runs were performed across each order.Data are reported as mean ± standard deviation.

Table 1
Computational times across PCE orders for the idealized artery model.Five PCE runs were performed across each order.Data are reported as mean ± standard deviation.